# Replication code for "Explaining the Oil Advantage" *World Politics* 67(2)
# Author: Paasha Mahdavi (paasha.mahdavi@gmail.com)
# 
# Data file required: iranreplicationdata.RData (or *.dta or *.csv)
# See "Data codebook.pdf" for details about variables and data sources
#
#
###############################################################################
#    Loading packages and data                                                #
###############################################################################

library(foreign); library(arm); library(texreg); library(ggplot2)
library(stargazer); library(nlme); library(reshape2); library(msm) 
library(plyr); library(xtable)

load("iranreplicationdata.RData")



###############################################################################
#                                                                             #
#-----------------------------------------------------------------------------#
#    TABLES                                                                   #
#-----------------------------------------------------------------------------#
#                                                                             # 
###############################################################################

###############################################################################
#    Summary Statistcs                                                        #
#    Table B.3                                                                #
###############################################################################

summarystats <- function(dbs){
  c(mean(dbs,na.rm=TRUE),
    sd(dbs,na.rm=TRUE),
    quantile(dbs,na.rm=TRUE,c(0,.25,.5,.75,1)),
    count(is.na(dbs))[2,2])
}

sss <- data[,c("incumbent","mining","lmining","sdistrict","prior","cleric","lngdpc.y","unemp","ethnic","public_emp","dpublic_emp_pct","bedsc","dbeds","stratio","dstratio")]

xtable((t(apply(sss,2,summarystats))),digits=c(0,rep(2,7),0))



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Table 2                                                                  #
###############################################################################

# Bivariate with time and interaction (column 1)
incRI.1 <- lme(incumbent ~  lmining*sdistrict + session, data, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding prior terms (column 2)
incRI.2 <- lme(incumbent ~  lmining*sdistrict + session + prior, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding Cleric dummy (column 3)
incRI.3 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric, data, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding logged GDP per capita (column 4)
incRI.4 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric + lngdpc.y, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding unemployment (column 5)
incRI.5 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric + lngdpc.y + unemp, data, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding ethnic minority dummy (column 6)
incRI.6 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric + lngdpc.y + unemp + ethnic, data, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Splitting sample into single member districts and multi member districts
# SMD (column 7)
incRI.single <- lme(incumbent ~  lmining + session + prior + cleric + lngdpc.y + unemp + ethnic, data[data$sdistrict==1,],na.action=na.omit,method="REML",  random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# MMD (column 8)
incRI.multi <- lme(incumbent ~  lmining + session + prior + cleric + lngdpc.y + unemp + ethnic, data[data$sdistrict!=1,], na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX table
vars <- c("Intercept","Resources (log)","SMD Dummy","Session (time)","Resources (log) $\times$ SMD","Number of prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic minority dummy")
texreg(list(incRI.1, incRI.2, incRI.3, incRI.4, incRI.5, incRI.6, incRI.single, incRI.multi), stars=c(0.01,0.05,0.10),reorder.coef=c(1,2,3,5,4,6,7,8,9,10),custom.coef.names=vars,digits=3)



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Robustness check #1: Using oil income pc instead of mining_pct           #  
#    Table B.4                                                                #
###############################################################################

inc.rob.1 <- lme(incumbent ~  lnoilpc*sdistrict + session, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

#RI models with covariates
inc.rob.2 <- lme(incumbent ~  lnoilpc*sdistrict + session + prior, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.3 <- lme(incumbent ~  lnoilpc*sdistrict + session + prior + cleric, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.4 <- lme(incumbent ~  lnoilpc*sdistrict + session + prior + cleric + lngdpc.y, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.5 <- lme(incumbent ~  lnoilpc*sdistrict + session + prior + cleric + lngdpc.y + unemp, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.6 <- lme(incumbent ~  lnoilpc*sdistrict + session + prior + cleric + lngdpc.y + unemp + ethnic, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
vars <- c("Intercept","Oil","SMD Dummy","Session (time)","Oil $\times$ SMD","Prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic dummy")
texreg(list(inc.rob.1, inc.rob.2, inc.rob.3, inc.rob.4, inc.rob.5, inc.rob.6), stars=c(0.01,0.05,0.10),reorder.coef=c(1,2,3,5,4,6,7,8,9,10),custom.coef.names=vars,digits=3)



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Robustness check #2: Using oil share instead of mining_pct               #   
#    Table B.5                                                                #
###############################################################################

#RI bivariate
inc.rob.1 <- lme(incumbent ~  lnoilshare*sdistrict + session, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

#RI models with covariates
inc.rob.2 <- lme(incumbent ~  lnoilshare*sdistrict + session + prior, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.3 <- lme(incumbent ~  lnoilshare*sdistrict + session + prior + cleric, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.4 <- lme(incumbent ~  lnoilshare*sdistrict + session + prior + cleric + lngdpc.y, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.5 <- lme(incumbent ~  lnoilshare*sdistrict + session + prior + cleric + lngdpc.y + unemp, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

inc.rob.6 <- lme(incumbent ~  lnoilshare*sdistrict + session + prior + cleric + lngdpc.y + unemp + ethnic, data,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
vars <- c("Intercept","Oil share","SMD Dummy","Session (time)","Oil share $\times$ SMD","Prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic dummy")
texreg(list(inc.rob.1, inc.rob.2, inc.rob.3, inc.rob.4, inc.rob.5, inc.rob.6), stars=c(0.01,0.05,0.10),reorder.coef=c(1,2,3,5,4,6,7,8,9,10),custom.coef.names=vars,digits=3)



###############################################################################
#    OLS  Modeling for Incumbency as DV                                       #
#    Robustness check #3: OLS with province and year fixed effects            #
#    Table B.6                                                                #
###############################################################################

data.na <- na.omit(data[,c("incumbent","lmining","sdistrict","session","prior","cleric","unemp","ethnic","province","lngdpc.y")])

# Bivariate with time and interaction (column 1)
inc.rob2.1 <- lm(incumbent ~  lmining*sdistrict + factor(session) - 1 + factor(province),data.na)

# Adding prior terms (column 2)
inc.rob2.2 <- lm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + factor(province),data.na)

# Adding Cleric dummy (column 3)
inc.rob2.3 <- lm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + factor(province),data.na)

# Adding logged GDP per capita (column 4)
inc.rob2.4 <- lm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + lngdpc.y + factor(province),data.na)

# Adding unemployment (column 5)
inc.rob2.5 <- lm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + lngdpc.y + unemp + factor(province),data.na)

# Adding ethnic minority dummy (column 6)
inc.rob2.6 <- lm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + lngdpc.y + unemp + ethnic + factor(province),data.na)

# With clustered standard errors (by province)
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) 
}

cl(data.na, inc.rob2.1, data.na$province)
cl(data.na, inc.rob2.2, data.na$province)
cl(data.na, inc.rob2.3, data.na$province)
cl(data.na, inc.rob2.4, data.na$province)
cl(data.na, inc.rob2.5, data.na$province)
cl(data.na, inc.rob2.6, data.na$province)

# LaTeX table
texreg(list(cl(data.na, inc.rob2.1, data.na$province), cl(data.na, inc.rob2.2, data.na$province), cl(data.na, inc.rob2.3, data.na$province), cl(data.na, inc.rob2.4, data.na$province), cl(data.na, inc.rob2.5, data.na$province), cl(data.na, inc.rob2.6, data.na$province)), stars=c(0.01,0.05,0.10), digits=3, omit="factor")

# LaTeX table without clustered standard errors
texreg(list(inc.rob2.1,inc.rob2.2,inc.rob2.3,inc.rob2.4,inc.rob2.5,inc.rob2.6),digits=2)



###############################################################################
#    Logit  Modeling for Incumbency as DV                                     #
#    Robustness check #4: Logit with province and year fixed effects          #
#    Table B.7                                                                #
###############################################################################

# Bivariate with time and interaction (column 1)
inc.rob3.1 <- glm(incumbent ~  lmining*sdistrict + factor(session) - 1 + factor(province),data.na, family=binomial(link="logit"))

# Adding prior terms (column 2)
inc.rob3.2 <- glm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + factor(province),data.na, family=binomial(link="logit"))

# Adding Cleric dummy (column 3)
inc.rob3.3 <- glm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + factor(province),data.na, family=binomial(link="logit"))

# Adding logged GDP per capita (column 4)
inc.rob3.4 <- glm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + lngdpc.y + factor(province),data.na, family=binomial(link="logit"))

# Adding unemployment (column 5)
inc.rob3.5 <- glm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + lngdpc.y + unemp + factor(province),data.na, family=binomial(link="logit"))

# Adding ethnic minority dummy (column 6)
inc.rob3.6 <- glm(incumbent ~  lmining*sdistrict + factor(session) - 1 + prior + cleric + lngdpc.y + unemp + ethnic + factor(province),data.na, family=binomial(link="logit"))

# LaTeX table
vars <- c("Intercept","Resources (log)","SMD Dummy","Resources (log) $\times$ SMD","Number of prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic minority dummy")
texreg(list(inc.rob3.1, inc.rob3.2, inc.rob3.3, inc.rob3.4, inc.rob3.5, inc.rob3.6), stars=c(0.01,0.05,0.10), omit.coef="factor",digits=3)



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Robustness check #5: REML with prov fixed effects & time random effects  #
#    Table B.8                                                                #
###############################################################################

# Bivariate with time and interaction (column 1)
inc.rob4.1 <- lme(incumbent ~  lmining*sdistrict + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding prior terms (column 2)
inc.rob4.2 <- lme(incumbent ~  lmining*sdistrict + prior + factor(province) - 1, data,na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding Cleric dummy (column 3)
inc.rob4.3 <- lme(incumbent ~  lmining*sdistrict + prior + cleric + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding logged GDP per capita (column 4)
inc.rob4.4 <- lme(incumbent ~  lmining*sdistrict + prior + cleric + lngdpc.y + factor(province) - 1, data,na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding unemployment (column 5)
inc.rob4.5 <- lme(incumbent ~  lmining*sdistrict + prior + cleric + lngdpc.y + unemp + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding ethnic minority dummy (column 6)
# This model can't converge because of singularity 
# between the ethnic minority dummy, time RE, and province FE
#
#inc.rob4.6 <- lme(incumbent ~  lmining*sdistrict + prior + cleric + lngdpc.y + unemp + ethnic + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX table
texreg(list(inc.rob4.1, inc.rob4.2, inc.rob4.3, inc.rob4.4, inc.rob4.5), omit="factor",digits=3)



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Robustness check #6: Removing resource-rich provinces                    #
#    Table B.9 (and B.10)                                                     #
###############################################################################

# Removing Khuzestan and Bushehr
data.rob <- data[data$province!="Khuzestan"&data$province!="Bushehr",]

# To run same analysis but with Kohgiluyeh and Boyerahmad removed
#data.rob <- data[data$province!="Bushehr"&data$province!="Kohgiluyeh and Boyerahmad",]
#data.rob <- data[data$province!="Khuzestan"&data$province!="Kohgiluyeh and Boyerahmad",]

# Bivariate with time and interaction (column 1)
inc.rob5.1 <- lme(incumbent ~  lmining*sdistrict + session, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding prior terms (column 2)
inc.rob5.2 <- lme(incumbent ~  lmining*sdistrict + session + prior, data.rob, na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding Cleric dummy (column 3)
inc.rob5.3 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding logged GDP per capita (column 4)
inc.rob5.4 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric + lngdpc.y, data.rob,na.action=na.omit,method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding unemployment instead of GDP (column 5)
inc.rob5.5 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric + lngdpc.y + unemp, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding ethnic minority dummy (column 6)
inc.rob5.6 <- lme(incumbent ~  lmining*sdistrict + session + prior + cleric + lngdpc.y + unemp + ethnic, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX table
vars <- c("Intercept","Resources (log)","SMD Dummy","Session (time)","Resources (log) $\times$ SMD","Number of prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic minority dummy")
texreg(list(inc.rob5.1, inc.rob5.2, inc.rob5.3, inc.rob5.4, inc.rob5.5, inc.rob5.6), stars=c(0.01,0.05,0.10),reorder.coef=c(1,2,3,5,4,6,7,8,9,10),custom.coef.names=vars,digits=3)



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Robustness check #7: Splitting into k-member districts                   #
#    Table B.11                                                               #
###############################################################################

mmd.1 <- lme(incumbent ~  lmining + session + prior + cleric + lngdpc.y + unemp + ethnic, data[data$seatsperdistrict==1,], na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE)) #single-member districts

mmd.2 <- lme(incumbent ~  lmining + session + prior + cleric + lngdpc.y + unemp + ethnic, data[data$seatsperdistrict==2,], na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

mmd.3 <- lme(incumbent ~  lmining + session + prior + cleric + lngdpc.y + unemp + ethnic, data[data$seatsperdistrict==3,], na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

mmd.4p <- lme(incumbent ~  lmining + session + prior + cleric + lngdpc.y + unemp + ethnic, data[data$seatsperdistrict>=4,], na.action=na.omit, method="REML", random= ~ 1 | province/district, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
vars <- c("Intercept","Resources (log)","Session (time)","Number of prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic minority dummy")
texreg(list(mmd.1,mmd.2,mmd.3,mmd.4p),stars=c(0.1,0.05,0.01),custom.coef.names=vars,digits=3)



###############################################################################
#    REML Modeling for Incumbency as DV                                       #
#    Robustness check #8: Data only for 2008 election                         #
#    Table B.12                                                               #
###############################################################################

# Data only for 2008
data.rob <- data[data$session==8,]

# Bivariate with time and interaction (column 1)
inc.rob6.1 <- lme(incumbent ~  lmining*sdistrict, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding prior terms (column 2)
inc.rob6.2 <- lme(incumbent ~  lmining*sdistrict  + prior, data.rob, na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding Cleric dummy (column 3)
inc.rob6.3 <- lme(incumbent ~  lmining*sdistrict  + prior + cleric, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding logged GDP per capita (column 4)
inc.rob6.4 <- lme(incumbent ~  lmining*sdistrict  + prior + cleric + lngdpc.y, data.rob,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding unemployment instead of GDP (column 5)
inc.rob6.5 <- lme(incumbent ~  lmining*sdistrict  + prior + cleric + lngdpc.y + unemp, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Adding ethnic minority dummy (column 6)
inc.rob6.6 <- lme(incumbent ~  lmining*sdistrict  + prior + cleric + lngdpc.y + unemp + ethnic, data.rob, na.action=na.omit, method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX table
vars <- c("Intercept","Resources (log)","SMD Dummy","Resources (log) $\times$ SMD","Number of prior terms","Cleric dummy","GDP per capita (log)","Unemployment rate","Ethnic minority dummy")
texreg(list(inc.rob6.1, inc.rob6.2, inc.rob6.3, inc.rob6.4, inc.rob6.5, inc.rob6.6), stars=c(0.01,0.05,0.10),custom.coef.names=vars,digits=3)



###############################################################################
#    REML Modeling for public goods as DV                                     #
#    All models with province fixed effects, time random effects              #
#    Table 3                                                                  #
###############################################################################

# DV: Public employment levels
patronageRI <- lme(public_emp*100 ~ lmining + lngdpc + factor(province) - 1, data, na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Public employment four-year changes
dpatronageRI <- lme(dpublic_emp_pct*100 ~ lmining + lngdpc + pub_emp + factor(province) - 1, data, na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Hospital beds levels
healthRI <- lme(bedsc ~ lmining + lngdpc + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Hospital beds four-year changes
dhealthRI <- lme(dbeds ~ lmining + lngdpc + bedstotal + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Student-teacher ratio levels
stratioRI <- lme(stratio ~ lmining + lngdpc + factor(province) - 1, data,na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Student-teacher ratio four-year changes
dstratioRI <- lme(dstratio ~ lmining + lngdpc + stratio + factor(province) - 1, data,na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Incumbent reelection (with no resources independent variable)
# Using lngdpc.y instead of lngdpc to include 1996 provincial GDP levels
inc.nooil <- lme(incumbent ~ session + prior + sdistrict + lngdpc.y + log(population) + public_emp + bedsc + stratio, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
texreg(list(patronageRI, dpatronageRI, healthRI, dhealthRI, stratioRI, dstratioRI, inc.nooil), stars=c(0.01,0.05,0.10), omit.coef="factor", digits=3)



###############################################################################
#    REML Modeling for public goods as DV                                     #
#    Robustness check #1: REML with province random intercepts                #
#    Table B.13                                                               #
###############################################################################

# DV: Public employment levels
patronage.rob1 <- lme(public_emp*100 ~ session + lmining + lngdpc, data, na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Public employment four-year changes
dpatronage.rob1 <- lme(dpublic_emp_pct*100 ~ session + lmining + lngdpc + pub_emp, data, na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Hospital beds levels
health.rob1 <- lme(bedsc ~  session + lmining + lngdpc, data, na.action=na.omit, method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Hospital beds four-year changes
dhealth.rob1 <- lme(dbeds ~  session + lmining + lngdpc + bedstotal, data, na.action=na.omit, method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Student-teacher ratio levels
stratio.rob1 <- lme(stratio ~  session + lmining + lngdpc, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Student-teacher ratio four-year changes
dstratio.rob1 <- lme(dstratio ~  session + lmining + lngdpc + stratio, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
texreg(list(patronage.rob1, dpatronage.rob1, health.rob1, dhealth.rob1, stratio.rob1, dstratio.rob1), stars=c(0.01,0.05,0.10), omit.coef="factor", digits=3)



###############################################################################
#    OLS Modeling for public goods as DV                                      #
#    Robustness check #2: OLS with prov fixed effects and year fixed effects  #
#    Table B.14                                                               #
###############################################################################

# DV: Public employment levels
patronage.rob2 <- lm(public_emp*100 ~ factor(session) + lmining + lngdpc + factor(province)-1, data, na.action=na.omit)

# DV: Public employment four-year changes
dpatronage.rob2 <- lm(dpublic_emp_pct*100 ~ factor(session) + lmining + lngdpc + pub_emp + factor(province)-1, data, na.action=na.omit)

# DV: Hospital beds levels
health.rob2 <- lm(bedsc ~  factor(session) + lmining + lngdpc + factor(province)-1, data, na.action=na.omit)

# DV: Hospital beds four-year changes
dhealth.rob2 <- lm(dbeds ~  factor(session) + lmining + lngdpc + bedstotal + factor(province)-1, data, na.action=na.omit)

# DV: Student-teacher ratio levels
stratio.rob2 <- lm(stratio ~  factor(session) + lmining + lngdpc + factor(province)-1, data,na.action=na.omit)

# DV: Student-teacher ratio four-year changes
dstratio.rob2 <- lm(dstratio ~  factor(session) + lmining + lngdpc + stratio + factor(province)-1, data,na.action=na.omit)

# LaTeX Table
texreg(list(patronage.rob2, dpatronage.rob2, health.rob2, dhealth.rob2, stratio.rob2, dstratio.rob2), stars=c(0.01,0.05,0.10), omit.coef="factor", digits=3)



###############################################################################
#    OLS Modeling for Public goods as DV                                      #
#    Robustness check #3: OLS with province fixed effects                     #
#    Table B.15                                                               #
###############################################################################

# DV: Public employment levels
patronage.rob3 <- lm(public_emp*100 ~ session + lmining + lngdpc + factor(province)-1, data, na.action=na.omit)

# DV: Public employment four-year changes
dpatronage.rob3 <- lm(dpublic_emp_pct*100 ~ session + lmining + lngdpc + pub_emp + factor(province)-1, data, na.action=na.omit)

# DV: Hospital beds levels
health.rob3 <- lm(bedsc ~  session + lmining + lngdpc + factor(province)-1, data, na.action=na.omit)

# DV: Hospital beds four-year changes
dhealth.rob3 <- lm(dbeds ~  session + lmining + lngdpc + bedstotal + factor(province)-1, data, na.action=na.omit)

# DV: Student-teacher ratio levels
stratio.rob3 <- lm(stratio ~  session + lmining + lngdpc + factor(province)-1, data,na.action=na.omit)

# DV: Student-teacher ratio four-year changes
dstratio.rob3 <- lm(dstratio ~  session + lmining + lngdpc + stratio + factor(province)-1, data,na.action=na.omit)

# LaTeX Table
texreg(list(patronage.rob3, dpatronage.rob3, health.rob3, dhealth.rob3, stratio.rob3, dstratio.rob3), stars=c(0.01,0.05,0.10), omit.coef="factor", digits=3)



###############################################################################
#    REML Modeling for public goods as DV                                     #
#    Robustness check #4: REML with province fixed effects and time random    # 
#                         effects and "singleprov" as control (see text)      #
#    Table B.16                                                               # 
###############################################################################

# DV: Public employment levels
patronage.rob4 <- lme(public_emp*100 ~ lmining + lngdpc + singleprov + factor(province) - 1, data, na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Public employment four-year changes
dpatronage.rob4 <- lme(dpublic_emp_pct*100 ~ lmining + lngdpc + singleprov + pub_emp + factor(province) - 1, data, na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Hospital beds levels
health.rob4 <- lme(bedsc ~ lmining + lngdpc + singleprov + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Hospital beds four-year changes
dhealth.rob4 <- lme(dbeds ~ lmining + lngdpc + singleprov + bedstotal + factor(province) - 1, data, na.action=na.omit, method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Student-teacher ratio levels
stratio.rob4 <- lme(stratio ~ lmining + lngdpc + singleprov + factor(province) - 1, data,na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# DV: Student-teacher ratio four-year changes
dstratio.rob4 <- lme(dstratio ~ lmining + lngdpc + singleprov + stratio + factor(province) - 1, data,na.action=na.omit,method="REML", random= ~ 1 | time, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
texreg(list(patronage.rob4, dpatronage.rob4, health.rob4, dhealth.rob4, stratio.rob4, dstratio.rob4), stars=c(0.01,0.05,0.10), omit.coef="factor", digits=3)



###############################################################################
#    REML Modeling with public goods as DV                                    #
#    Extra incumbency models for Table 3, Column 7                            #
#    (Not in article or appendix)                                             #
###############################################################################

# Health and education levels
inc.nooil.1 <- lme(incumbent ~ session + prior + sdistrict + bedsc + stratio, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Health and education changes
inc.nooil.2 <- lme(incumbent ~ session + prior + sdistrict + dbeds + dstratio, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Public employment levels
inc.nooil.3 <- lme(incumbent ~ session + prior + sdistrict + public_emp, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Public employment changes
inc.nooil.4 <- lme(incumbent ~ session + prior + sdistrict + dpublic_emp, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# PG levels
inc.nooil.5 <- lme(incumbent ~ session + prior + sdistrict + public_emp + bedsc + stratio, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# PG changes
inc.nooil.6 <- lme(incumbent ~ session + prior + sdistrict + dpublic_emp + dbeds + dstratio, data,na.action=na.omit,method="REML", random= ~ 1 | province, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# LaTeX Table
texreg(list(inc.nooil.1,inc.nooil.2, inc.nooil.3, inc.nooil.4, inc.nooil.5, inc.nooil.6),stars=c(.1,.05,.01), digits=4)



###############################################################################
#    Checking correlation with resources and turnout                          #
#    Using province-level data                                                #
#    Table B.18                                                               #
###############################################################################

# Loading and merging turnout data
# NOTE: requires provincedata.csv and turnout.csv
provincedata <- merge(read.csv("provincedata.csv"),read.csv("turnout.csv"),by=c("province","session"),all.x=TRUE)

# Visualizing correlation with resources
with(provincedata,plot(lmining,turnout))

# Computing Pearson correlation with resources
with(provincedata,cor(lmining,turnout,"pairwise.complete.obs"))
# 0.274

# Computing Pearson correlation with incumbent reelection
with(provincedata,cor(incumbent,turnout,"pairwise.complete.obs"))
# 0.108

# OLS coefficients for resources on turnout
summary(lm(turnout~lmining+factor(province)-1,provincedata))
# Beta: 2.7, std err: 2.8, p-value: 0.34

# OLS coefficients for incumbent reelection on turnout
summary(lm(turnout~incumbent+factor(province)-1,provincedata))
# Beta: 8.7, std err: 9.5, p-value: 0.37



###############################################################################
#    Hausman test                                                             #
#    For province random effects vs fixed effects                             #
#    Footnote 74 in manuscript                                                #
###############################################################################

# For simplest model with year dummies
require(plm)
form <- incumbent ~ lmining + sdistrict + lmining*sdistrict + factor(session) - 1
wi <- plm(form, data = data.na, index = c("province"), model = "within")
re <- plm(form, data = data.na, index = c("province"), model = "random")
phtest(wi, re)





###############################################################################
#                                                                             #
#-----------------------------------------------------------------------------#
#    FIGURES                                                                  #
#-----------------------------------------------------------------------------#
#                                                                             # 
###############################################################################

###############################################################################
#    Candidate Vetting                                                        #    
#    Figure 1                                                                 #
###############################################################################

vet <- data.frame(cbind(year=seq(1984,2008,by=4),candidates=c(1584,2001,3150,5365,6849,8145,7597),eligible=c(1161,1615,2090,3276,6083,4561,4476),incumbents=c(0,0,40,NA,NA,NA,48)))

p <- ggplot(vet)
p + geom_rect(aes(xmin=year-1,xmax=year,ymin=0,ymax=candidates)) + geom_rect(aes(xmin=year,xmax=year+1,ymin=0,ymax=eligible),fill="grey") + geom_point(aes(x=year-.5, y=candidates)) + geom_line(aes(x=year-.5, y=candidates)) + geom_point(aes(x=year+.5, y=eligible),col="grey") + geom_line(aes(x=year+.5, y=eligible),col="grey") + theme_bw() + xlab("Election Year") +  ylab("Number of Candidates") + scale_x_continuous(breaks=seq(1984,2012,4))



###############################################################################
#    Map of resources (done in QGIS)                                          #
#    Figure 2                                                                 # 
###############################################################################

# For the GIS map of average mining pct across all years, done in QGIS
formap <- log(100*tapply(data$mining,list(data$province),mean, na.rm=TRUE))
write.csv(formap,'formap.csv')



###############################################################################
#    Predicted interaction of SMD dummy and resources variable                #
#    Figure 3                                                                 #
###############################################################################

# main model but using lmer instead of lme
fit1 <- lmer(incumbent~session + lmining*sdistrict + prior + cleric + lngdpc.y + unemp + ethnic+ (1|district/province), data=data, na.action=na.omit)

# NA-omitted data for constructing predicted plot
incoildata <- na.omit(data[,c("incumbent","session","lmining","sdistrict","prior","cleric","lngdpc.y","unemp","ethnic","province","district")])
incoildata$lmining.sdistrict <- incoildata$lmining*incoildata$sdistrict
incoildata$intercept <- 1
mu.x <- as.matrix(apply(incoildata[,c(13,2:9,12)],2,mean))
session <- rep(mu.x[2],22)
lmining <- rep(seq(-10,0,by=1),2)
sdistrict <- rep(0:1,1,each=11)
prior <- rep(mu.x[5],22)
cleric <- rep(mu.x[6],22)
lngdpc <- rep(mu.x[7],22)
unemp <- rep(mu.x[8],22)
ethnic <- rep(mu.x[9],22)
newdata1 <- data.frame(intercept=rep(1,22),session=session,lmining=lmining,sdistrict=sdistrict,prior=prior,cleric=cleric,lngdpc=lngdpc,unemp=unemp,ethnic=ethnic,lmining.sdistrict=lmining*sdistrict)
newdata.m <- as.matrix(newdata1,ncol=10,byrow=T)
estmeans <- fit1@beta
newdata1$mineP <- newdata.m%*%matrix(estmeans)

# Figure 3
pdf(file="predictedinteraction.pdf",width=10,height=8)
plot(newdata1$lmining[1:11],newdata1$mineP[1:11],axes=FALSE,xlim=c(-10,0),ylim=c(0,.6),xlab="Oil and minerals value-added (% Resource GDP)",ylab="Predicted probability of incumbent reelection",type='n',cex.lab=1.5)
lines(newdata1$lmining[12:22],newdata1$mineP[12:22],lwd=3)
lines(newdata1$lmining[1:11],newdata1$mineP[1:11],lwd=3,lty=2)
points(jitter(data$lmining,1), jitter(rep(0.05,length(data$lmining)),10), pch=20, cex=.5)
axis(1,at=lmining,labels=round(100*exp(lmining),2),cex.axis=1.5)
axis(2,at=seq(.0,.6,.1),labels=seq(0,60,10),cex.axis=1.5)
text(x=-3,y=.5,"SMD",cex=1.5)
text(x=-3,y=.2,"MMD",cex=1.5)
dev.off()



###############################################################################
#    Partial plot for vote share regression                                   #
#    Figure 4                                                                 #
###############################################################################

votedata <- na.omit(data[,c("vote","incumbent","lmining","unemp","lngdpc.y","cleric","session","prior","province","sdistrict")])

# main model for vote shares
vote.1 <- lm(vote ~ lmining + unemp + lngdpc.y + cleric + session + prior + factor(province)-1,data=votedata, subset=incumbent==1&sdistrict==1)
summary(vote.1)

# residuals
vote.resid <- lm(vote ~ unemp + lngdpc.y + cleric + session + prior + factor(province) - 1, data=votedata, subset=incumbent==1&sdistrict==1)
oilvote.resid <- lm(lmining ~ unemp + lngdpc.y + cleric + session + prior + factor(province) - 1, data=votedata, subset=incumbent==1&sdistrict==1)
vote.main.resid <- vote.resid$residual
oilvote.main.resid <- oilvote.resid$residual

# Figure 4
pdf(file="voteshare.pdf",width=10,height=8)
p <- data.frame(oilvote.main.resid, vote.main.resid)
ggplot(aes(x=oilvote.main.resid,y=vote.main.resid),data=p) + geom_point(size=3) + geom_smooth(method="lm",size=1,col="black") + theme_bw() + xlab("Oil & Minerals Value Added, logged | controls") +  ylab("Incumbent vote share | controls") + theme(axis.title = element_text(size = 22, vjust = .5), axis.text = element_text(size=20))
dev.off()



###############################################################################
#    Partial plot for student-teacher ratio levels and changes regressions    #
#    Figure 5                                                                 #
###############################################################################

datastratio <- data[,c("province","session","stratio","lmining","lngdpc")]
datastratio <- na.omit(datastratio)

# main model for student-ratio levels
stratioRI.main <- lme(stratio ~ lngdpc + lmining + factor(province) - 1, datastratio, na.action=na.omit, method="REML", random= ~ 1 | session, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# residuals
stratioRI <- lme(stratio ~  lngdpc + factor(province) - 1, datastratio,na.action=na.omit,method="REML", random= ~ 1 | session, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))
stratio.resid <- lme(lmining ~ lngdpc + factor(province) - 1, datastratio,na.action=na.omit,method="REML", random= ~ 1 | session,control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

str.resid <- stratioRI$residual[,1]
oil.resid.str <- stratio.resid$residual[,1]

# Figure 5 panel a
pdf(file="stratio_partial.pdf",width=10,height=8)
p <- ggplot(aes(x=oil.resid.str,y=str.resid),data=datastratio)
p + geom_point(size=3) + geom_smooth(method="lm",size=1,col=1) + theme_bw() + xlab("Oil & Minerals Value Added, logged | controls") +  ylab("Student-Teacher Ratio | controls") + theme(axis.title = element_text(size = 22, vjust = .3), axis.text = element_text(size=20)) + scale_y_continuous(breaks=seq(-3,6,1.5))
dev.off()

# -- Now looking at changes in student-teacher ratio using dstratio --------- #

datadstratio <- data[,c("province","session","stratio","dstratio","lmining","lngdpc")]
datadstratio <- na.omit(datadstratio)

# main model for student-teacher ratio changes
dstratioRI.main <- lme(dstratio ~ lngdpc + lmining + stratio + factor(province) - 1, datadstratio, na.action=na.omit, method="REML", random= ~ 1 | session, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# residuals
dstratioRI <- lme(dstratio ~  lngdpc + stratio + factor(province) - 1, datadstratio,na.action=na.omit,method="REML", random= ~ 1 | session, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))
dstratio.resid <- lme(lmining ~ lngdpc + stratio + factor(province) - 1, datadstratio,na.action=na.omit,method="REML", random= ~ 1 | session,control=lmeControl(msMaxIter = 200, msVerbose = TRUE))
dstr.resid <- dstratioRI$residual[,2]
oil.resid.dstr <- dstratio.resid$residual[,2]

# Figure 5 panel b
pdf(file="dstratio_partial.pdf",width=10,height=8)
p <- ggplot(aes(x=oil.resid.dstr,y=dstr.resid),data=datadstratio)
p + geom_point(size=3) + geom_smooth(method="lm",size=1,col=1) + theme_bw() + xlab("Oil & Minerals Value Added, logged | controls") +  ylab("4-yr Changes in Student-Teacher Ratio | controls") + theme(axis.title = element_text(size = 22, vjust = .3), axis.text = element_text(size=20)) +   scale_y_continuous(breaks=seq(-.25,.25,.05))
dev.off()



###############################################################################